Mon. Not. R. Astron. Soc. 000. UMTS! (20041 Printed 2 February 2008 (MN WF&L style file v2.2) 



Time-dependent secular evolution in galaxies 

Martin D. Weinberg* 

Department of Astronomy, University of Massachusetts, Amherst, 
2 February 2008 

ABSTRACT 

iLvnden-Bell fc Kalnai d l)l972l LBK) presented a useful formula for computing the 
long-range torque between spiral arms and the disk at large. The derivation uses 
| second-order perturbation theory and assumes that the perturbation slowly grows 

. over a very long time (the time- asymptotic limit). This formula has been widely used 

^ ( ' to predict the angular momentum transport between spiral arms and stellar bars 

between disks and dark-matter halos. However, this paper shows that the LBK time- 
asymptotic limit is not appropriate because the characteristic evolution time for galax- 
ies is too close to the relevant dynamical times. We demonstrate that transients, not 
present in the time-asymptotic formula, can play a major role in the evolution for 
realistic astronomical time scales. A generalisation for arbitrary time dependence is 
£> ' presented and illustrated by the bar-halo and satellite halo interaction. The natural 

on : time dependence in bar-driven halo evolution causes quantitative differences in the 

overall torque and qualitative differences in the physical- and phase-space location of 
^-j- . angular momentum transfer. The time-dependent theory predicts that four principal 

f^S ' resonances dominate the torque at different times and accurately predicts the results 

of an N-body simulation. In addition, we show that the Inner Lindblad Resonance 
(ILR) is responsible for the peak angular momentum exchange but, due to the time 
. dependence, the changes occur over a broad range of energies, radii and frequencies. 

Qh We describe the implication of these findings for the satellite-halo interaction using a 

q , simple model and end with a discussion of possible impact on other aspects secular 

5— ( ' galaxy evolution. 
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1 INTRODUCTION 

A near-equilibrium galaxy only evolves through the excitation of some non-axisymmetric structure such as spiral arms, bars, 
infalling satellites, etc. The resulting "waves" are likely actors in triggering star formation, mediating inward gas flow thereby 
fueling AGN, heating the stellar disk, among other observable phenomena. In addition, the overall evolution of near equilibrium 
galaxies is caused by the collective transport of energy and angular momentum from these same structures. An accurate theory 
of secular evolution is essential for developing our understanding of these interrelated dynamics. 

The first step is predicting the evolution of the non-axisymmetric distortions. The underlying dynamical mechanism for 
the stellar and dark matter components is understood. At the level of individual orbits, an excitation such as a bar, a spiral 
pattern, or orbiting satellite produces a periodic distortion on all orbits, much like a driven harmonic oscillator. Although the 
torque distorts a typical orbit, this distortion averages to zero over many orbital periods for most orbits. However, for nearly 
closed orbits 1 , the net torque does not vanish. An adiabatic invariant is broken by the vanishing precession frequency of the 
closed orbit; clearly, there is no adiabatic invariant for a degree of freedom with no motion. The perturbation then couples 
to the oscillatory distortion for these orbits and net torque is transferred to or from these orbits. This coupling between the 
forced oscillation and the perturbation will be second order in the perturbation strength. In time, the perturbation itself and 

* E-mail: weinberg@astro.umass.edu 

1 A closed orbit is one with integral commensurabilities between the orbital frequencies and the pattern frequency. 
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the underlying equilibrium is slowly changed by this torque. This lea ds to a finite measu re of orbits that transfer angular 
momentum to and from the pattern (for additional discussion, see IWeinberg fc Katzll2004l. hereafter Paper 2). 

Although the first-order perturbation may dominate the instantaneous changes to phase space, the second-order changes 
describe the net changes tha t will persist after many dynamical times. This process was described mathematically by 
iLvnden-Bell fc Kalnaisl (|l972l hereafter LBK) in the limit that the change in the perturbation is very slow but still fast 
enough that the orbital perturbations remain linear. LBK derived a formula describing the exchange of angular momentum 
between a spiral pattern and the rest of the disk assuming that the perturbation began infinitely long in the past, the time- 
asymptotic limit. In this limit, any vestigial respon se from the formation of the pattern is gone. The LBK app r oach has been 
widel y applied to estimate s ecular evolution (e.g. iGoldreich fc Tremaineiri979l : IWeinberall985l : IZhandllQQSi : lAthanassoulal 
l2003h . ICarlberg fc Sellwoodl lll985l) took the first step towards a time-dependent generalisation of the LBK formalism for 
studying secular disk evolution. They included time-dependent transients but evaluated the secular changes in the infinite 
time limit. This is appropriate for studying the cumulative effect of short-lived transients. However, the perturbation theory 
may be fully generalised to treat arbitrary time-dependent perturbations over arbitrary intervals of time. This allows us to 
treat the long-lived bar and satellite interactions over the full age of a galaxy with some surprising results. 

In attempting a detailed description of the angular momentum exchange between a rotating bar and a dark halo in an N- 
body simulation (see Paper 2), I found a significant discrepancy between the predictions of the second-order LBK perturbation 
theory and the N-body simulation. This discrepancy persisted for perturbations over a wide range of bar amplitudes and scaled 
with time and perturbation amplitude as expected from second-order linear theory, implying that its source is not a break 
down in linearity. Rather, the problem is due to the assumption of an infinitely slow growth of the perturbation. Not only 
is the number of characteristic dynamical times in a galactic age modest, the growth of a bar, arms and of course satellites 
is most likely a small fraction of the galactic age. Because patterns often appear over several orbital periods and the total 
number of orbital periods available are small, transients may be significant and the time-asymptotic limit does not apply. I will 
show in this paper that a finite-time-limit generalisation of the LBK formula gives quantitatively and qualitatively different 
results. Quantitatively, we will see that the overall torque is smaller than that computed from the LBK formula for a rotating 
bar. Qualitatively, different commensurabilities are important at different epochs of the evolution. Because the location of 
the resonances in physical space are important for the long-term evolution of galaxy, including the finite-time response is 
important for understanding this evolution. In this limit, the time dependence is more than simply a transient but embodies 
the galaxy's evolutionary history. 

I will present a new, generalised secular evolution formula and describe the important details of the dynamics in J21 
The new torque formula replaces the delta function in the LBK formula by an integral over the time-dependent perturbation 
and is similarly straightforward to apply. Two specific applications are presented in bar-dark halo coupling and satellite- 
dark halo coupling. N-body simulations will demonstrate the importance of the finite-time response to the net torque and 
illustrate the discrepancy with the LBK formula. In H3.ll we will see that the ILR dominates the bar evolution and use this 
new formalism to explicitly predict the location of the angular momentum deposited. The same dynamics is then applied to 
sinking satellites in M.2I We will conclude with 21 



2 BASIC PRINCIPLES 

The total torque on an ensemble of orbits perturbed from equilibrium is 

where F is the phase-space domain and /(x, v, i) is the phase-space distribution function, L z = dL z /dt is the torque per orbit, 
and (•) denotes the phase-space average. The angular momentum L z is a conserved quantity in an equilibrium axisymmetric 
system and therefore L z is first-order in the perturbation amplitude. The distribution function has, of course, non-vanishing 
zeroth- and first-order terms. Over a long period of time, the first-order contribution to (L z ) averages to zero. At the next 
order, the oscillating components from the first-order dependence in (L z ) and the first-order induced change in the distribution 
function, /i(x, v, t), coherently reinforce each other, leading to a non- vanishing contribution as t — > oo. However, a closed-form 
solution to this problem does not demand the t — > oo limit. We will derive the first-order changes in these two quantities in 
H2.1l and H2.2l and put these together H2.3l to derive the torque for a finite duration perturbation. 

2.1 Change in angular momentum for a single orbit 

Hamilton's equations determine the rate of change in angular momentum for each orbit as follows: 
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where [■, ■] are the Poisson brackets. Because L z is a conserved quantity in the absence of any perturbation, this equation 
becomes 

dL z _ dH 9L ± _dH_ dL ± ^_<Mi dL z 2 

dt dl ' <9w dw' dl dw ' dl [ ' ( ' 

using action-angle variables I, w, Hi is the first-order, perturbed Hamiltonian and A is the perturbation amplitude. The 
subscripts '0' and 'f ' will indicate zeroth- and first-order quantities, respectively. 

The orbits in an axisymmetric equilibrium galaxy are purely quasi-periodic in the angles w and the actions I are invariant. 
The advantage of action-angle variables is that any phase-space quantity may be expanded in a Fourier series in w. In 
particular, we now expand Hi in action-angle variables: 

Hi(I,w,t)=^i?n(I J i)e il ' w (4) 
1 

where 

ff„(I,t)=^3| dw e - llw ffi(I,w,t) (5) 
and 1 = (h,h, h) is a vector of integers, one for every degree of freedom. Applying the angle transform to equation yields JjU 

^ = -5>3tfii(MK'' w (6) 
i 

Physically, equations © and @ describe the first-order evolution of a particular orbit's z angular momentum component. 



2.2 Changes in the distribution function 

The development begins with the first-order solution of the linearised collisionless Boltzmann equation (CBE): 
dfi dH df dH x df _ 

dt dl dw dw dl 1 ' 

The quantities /o and Hq are the unperturbed phase-space distribution function and the unperturbed Hamiltonian, etc. We 
will solve for the phase-space dependence by an action-angle expansion (as in eq. |1J and for the time evolution by Laplace 
transform of the CBE. We denote a Laplace transformed variable by a tilde and an action-angle transformed variable by 
subscript 1. The Fourier-Laplace transform of the CBE is 

s/i+il-fi/i-il-^tfi^O (8) 
or, upon solving for /i, 

h{l,s) = —f-^Hi^s). (9) 
s + ii ■ Si 

where ft = dH /dI. We will now show that the Laplace transform rather than inclusion of a term exp(ei) the limit e — > 
to impose the arrow of time allows one to solve for an arbitrary time dependence. The total perturbing potential Hi is 
the response of the galaxy combined with the external perturbation. The quantity flu (I, s) is the Laplace transform of the 
time-dependent coefficient of the action-angle expansion describing the perturbation: 



/•oc 

Hu(I,s)= / dte- st Hn(I,t). 
Jo 



(10) 



In cases where the perturbation is convenie ntly described in spherical harmonics, intrinsic symmetries allow Hi i to be further 
simplified fsee iTremaine fc Weinberslll984l) . 

The desired solution is the inverse Laplace transform of the series: 



oo 



where each coefficient /i is given by equation Q. The inverse Fourier-Laplace transform of equation il l H is straightforward 
albeit cumbersome. Recall that the product of the first-order distribution function and the first-order torque per orbit gives 
the second-order torque. To proceed we first need to transform equation © back to the time domain. Combining equations 
itTTt and 0, tne 

inverse Laplace transform of equation © requires evaluation of the following expression: 

1 pc+ica a t pea pea ^ pc+iaa s (i-i') 

J c — ioc J J J c — zoo 
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/ dt'Hu(I,t' 
Jo 

{ dt'Hn(I,t')e 
Jo 



e -a-nct-f) iit>t , 
otherwise 



i^jhn(t'-t) ^2) 

where c is more positive than the any exponent of divergence in Hn(I,t). This is not restrictive for any realistic astronomical 
perturbation. The second equality follows by exchanging the order of integration and using Cauchy's theorem. Putting this 
together, we have the following time-dependent solution for the first-order distribution function 



A(w,I,t)= il.|e fl -|r ( it'i/ 11 (I,t')e M ' t, - t »] 

7, 7., J„.— WO ) 



(13) 



2.3 Second-order torque 

To get the second-order torque as a function of actions (or other conserved quantities following Jeans' theorem), we multiply 
by equation © by equation I13H and integrate over all of phase space: 

(*h\ = fdif dw [/ (I) +/l (I, w) + •■■]= [dlf dw^/ 1 (I,w) + 0(A 3 ). 



\ dt I J J dt ' ° ' J J dt 

The first term is oscillatory in w and vanishes. The product of the two first-order solutions, denoted here by subscripts 1 and 
1', has an angle dependence of the form 

exp(— il ■ w) exp(il ■ w) = exp(i[I — l'] ■ w) 

and therefore the angle integral is only non-zero when 1' = 1. Integrating this product over the remaining phase-space variables 
(actions) gives an expression for the total torque on the phase space from the perturbation: 



(2tt) 3 J2 ls f t dt'H 11 (I,t')e il ^ t '-AH 1 - l (I,t). 



(14) 



( dL, 
\ dt 

H ,13= — 00 

Since physical quantities are real, equation @ implies H*i(l, t) = Hi-\(1, t). We may now change variables from actions to 
energy and angular momentum and the third action L z , the 2-component of the angular momentum, in terms of its ratio to 
the total angular momentum J as follows: cos/3 = L z /J. The angle (3 describes the colatitude of the orbital plane. Equation 
1141 1 becomes 

For a fixed gravitational potential profile whose only time dependence is through a rotating pattern, the time dependence 
in the perturbation separates as _ffi(I,w,t) = a(t)Hi(I,w) and therefore Hn(I,t) = a(t)Hn(T). Furthermore for a constant 
pattern speed, the time-dependent coefficients have the form: a(t) = a(0) exp(imf2 p i). The time-dependent terms in equation 
115H then take the simple form 

dt'e a - 0(t '- t) a(t')|o*(t) = J <ft'e a ' nt 'a(*')W(*) 

sin[(l ■ ft — mSl p )t] 



= a(Q)a*(0)- 



1 • f2 — mS7 p 



— > a(0)a*(0)7r5(l ■ Ct - mVL v ) (16) 

where the final expression obtains in the limit t —* 00. Upon substituting into equation 11511 . we recover the LBK formula. 
The mathematics of this convergence also indicates a possible numerical pitfall in applying equation 1141 : at large times T, 
equation 116^ will oscillate rapidly with E and J or I and may require special treatment. 

Equation l|15|l can be integrated to find the run of angular momentum L z with time for some arbitrary time-dependent 
perturbation. At each time step, one must expand Hi(I,t) for the desired range of 1, integrate over the temporal history, 
and then integrate over phase space. The solution for a general time dependence is computationally intensive, therefore. In 
the slowing-bar example ( A3. It . the problem has separable time dependence of the form Hi(I,w,t) = a(t)H\(l, w). In this 
equation 1151 may be written 

' dL z \ _ /r) _^3 , f f f dEd,JJd(cos f3) 1 df a 



iai^ z \ . .3 v-^ 1 dJjydJ j (lycos p ) 
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jjf dt'c a - n ( t '- t) o(t')|a*(t)fl-ii(I)flJ_ I (I). (17) 

In this case, a grid of Hii(I) along with the other phase-space dependent terms may be determined once and tabled for 
all time, simplifying the calculation. The integral J Q di'exp[il • Sl(t' — t)]Hn(l,t') must be computed with sufficient small 
time steps to resolve the oscillations obtained with variations in I. Finally, the right- and left-hand sides of equation 11411 are 
generally coupled. For a bar-halo interaction with a fixed bar shape for example, the patter n speed is O p = L z /I z where J z is 
the bar's moment of inertia. Stability of the solution is then a concern. I used Gear's method loeail lll96flF) to maintain stability 
for the numerical solutions presented here. The problem is numerically more stable when the right hand side is independent of 
L z as is the case for a satellite-halo interaction whose orbital decay is computed from local Chandrasekhar dynamical friction 



3 EXAMPLES 
3.1 Slowing bar 

Following Paper 2, we compare the resu lts of N-body simulations and linear theory for an ellipsoidal bar in an NFW dark- 
matter halo profile iNavarro et aljll997l) : p(r) oc r _1 (r + r 3 )~ 2 . The fiducial model has a bar of length (semi-major axis) r a 
with and axis ratios b/a = 0.2, c/a — 0.05 and bar mass equal to one-half of the enclosed dark-halo mass. The total mass of 
the bar is 0.1 of the halo mass inside the virial radius for a concentration c = 15. A realistic present-day stellar bar would have 
a length roughly 25% of r s . However, this length greatly increases the number of particles required in the N-body simulation. 
Because the NFW profile is scale free inside of r s , the behaviour should scale to smaller bars. This is verified in Paper 2. 

The pattern speed evolves while conserving the total angular momentum of the initial bar and dark halo system. The 
bar perturbation is slowly turned on according to f(t) = (1 + erf[(t — Tq)/A])/2 with To = 1/2 Gyr and A = 1/4, which 
corresponds to a full-width turn-on time of 1 Gyr scaled to the Milky Way. Bar formation in simulations occurs over a small 
number of dynamical timescales (e.g. within 200 Myr scaled to the Milky Way). An assumption of To — 1/2 and A = 1/4 
then is an upper bound to a realistic formation time scale; we are choosing the slowest time dependence possible, accounting 
for bar length scale. Paper 2 describes the additional details of the simulation procedure. 

Figure Q compares the resulting evolution of the N-body simulation, the linear theory presented in an d the LBK 
formula. Because the bar figure is constant in time, the torque takes the form of equation I17H . The N-body simulation 
closely follows the time-dependent linear theory prediction. The evolution is significantly slower than predicted by the LBK 
formula! Clearly, the time-dependence changes the t orque. The second panel of Figure repeats the simulation and linear 
theory prediction for a bar in a King model profile (jKin elll96fj) chosen to match the cumulative mass profile of the NFW 
profile inside of the virial radius (Wo — 6, rr — 2r V i r ) demonstrating that the differences due to time dependence are not a 
consequence of the NFW profile. 

Equation l|15|l or l|17|l allows us to determine the contributions from each resonance separately (see Fig. |2J and determine 
the origin of the differences between the time-dependent theory and LBK. The qualitative differences between the case with 
explicit time-dependence and the time-asymptotic LBK regime are evident. For the LBK case, the torque is dominated by 
the ILR-like 1 = ( — 1,2,2) resonance. The contribution peaks as the bar begins to slow and decreases when the pattern 
frequency becomes too small to easily couple to low-order resonances in the vicinity of the bar. For the time-dependent case, 
the corotation resonance (0,2,2) dominates at early times. During the rapid slowing phase, 0.4 < t < 1, the ILR ( — 1,2,2) 
dominates. The ILR is weaker in the time-dependent theory and the simulation. This accounts for a large fraction of the 
difference in Figure At later times, a near equilibrium obtains between the various commensurabilities, some without true 
resonances. These all converge to the LBK limit as t — ► oo. 

In Figure [3] I show the torque from the ILR summed over all angular momentum plotted as a function of energy E 
(top axis) and radius r of the circular orbit with the same energy E and time t (bottom axis) . The peak angular momentum 
deposition occurs inside the bar radius (due to the ILR) and corresponds with the radial location of the profile change 
described in Paper 2. Comparison of Figures and confirms that the peak torque occurs during the rapid pattern speed 
evolution which also accounts for the broad region of this contribution in energy for 0.4 < t < 1. For t > 1, the ILR oscillates 
at lower amplitude as it approaches the constant t — > oo LBK value. 

We now take a closer look at the contributions from all terms (h,h,l3 = m) in equation 11711 . The contribution to the 
torque for each term and the total is shown in Figure for a variety of different bar lengths and turn-on parameters To, A. 
The mass of the bar has been decreased by a factor of five to Mb = 0.01 to lengthen the slow down period and help isolate 
the contribution of each resonance. For clarity, only the top four contributing resonances are drawn. The first panel of the 
describes the torque for a bar with length r s with initial period 1/3, To = 1/2 and A = 1/4; the bar grows to full strength 
over several bar periods. At t < 1/2, early in the bar growth, the corotation (0, 2, 2) and the (1, 0, 2) resonances dominate the 
torque (see also Fig.[5Jl. At the peak of the pattern-speed decrease, t « 0.7, the ILR dominates the torque. Moreover, the ILR 
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(a) NFW 




(b) King 



Figure 1. The pattern speed evolution of a bar in a NFW (a) and King model (b) dark-matter halo for an N-body simulation and 
computed from the second-order perturbation theory in the time-dependent and LBK limits. The time-dependent theory accurately 
follows the patter speed evolution while the time-asymptotic LBK theory does not. 
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Figure 2. Relative contributions to the total torque L z = I z Cl z from each resonance for time-asymptotic (LBK) regime (a) and with 
explicit time-dependence (b). Each curve in the figure shows a single 1 term from equation 1151 as labelled. The run of il z derived from 
the total torque is also shown. 
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Figure 3. ILR contribution to the torque (Q z = L z /I z ) as a function of energy (top axis) and corresponding radius of the corresponding 
circular orbit (bottom axis) as a function of time as the pattern speed evolves. This figure is obtained by solving equation 1151 and 
plotting the energy contribution to the phase-space integral as a function of time for the same example shown in Fig. [5] 



dominates to the integrated torque for the entire slowdown event. Note that terms 1 without pattern-speed commensurabilities, 
iiOi + ^2^2 7^ mQp, through the phase-space distribution also contribute transiently but significantly to the total torque. For 
example, the (1,-2,2) term, an "negative" ILR contributes to the torque but would only be commensurate for a negative 
pattern speed. This transient echo decays after the rapid slow down phase since it is purely oscillating without a net second- 
order contribution. Nonetheless, it diminishes the torque during peak of the slow down. The ILR and the other commensurate 
terms, on the other hand, approaches the finite non-zero time-asymptotic value predicted by the LBK formula. 

The relative importance of the torque to the inner halo is nicely visualised by computing the net change in angular 
momentum inside of some energy E from equation l|150 as a function of time. That is, let L(< E) denote equation l|15|l with 
the energy integral from [E m in, E]. Then, define the cumulative change of angular momentum in time as 

AL(< E,<t) = J dtt(< E) (18) 

where t is the initial time. The total orbital angular momentum in the initial phase space with energy less than E is 

S« E) = (2nf HI ^^f m J) J- (19) 

Then, the ratio of these two quantities, AL(< E, < £)/£(< E) indicates the importance of the resonant torque to the structure 
of the inner halo profile. This ratio is shown in Figure^Jfor the bar-halo interaction depicted in Figure[3]for the ILR resonance. 
The contour scale shows the fractional change of the angular momentum in the halo inside of a given energy or radius. A 
large fractional change indicates a strong possibility of significant structural change. The peak absolute torque from the ILR 
occurs near r « 0.02 and t « 0.6 (Fig.0 while the peak relative torque occurs near r w 0.003 (Fig. 0. In particular, N-body 
simulations from Paper 2 show that this bar-halo interaction will flatten the NFW cusp for 0.001 < r < 0.01 as predicted by 
Figure |1] Hence, this comparison suggests that a 15%- 20% fractional change is sufficient to drive inner-halo evolution. 

The sequence in Panels (a)-(c) in Figure [£] shows the same bar model with different turn on parameters To and A. The 
corresponding distribution of torque deposited as a function of time and energy (radius) is similar to Figure [3] For longer 
intervals of bar growth, the resonances are more localised in energy and radius. The steep ILR peak during the rapid evolution 
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Figure 4. Ratio of the cumulative angular momentum contributed by the ILR inside of energy E and radius r of the corresponding 
circular orbit at time t to the total angular momentum of the halo inside of energy E. In other words, the contour scale shows the 
fractional change of the angular momentum in the halo from the ILR resonance inside of a given energy or radius. 

phase covers a broad swath in energy and radius and leads to profile evolution. At later times as the system approaches the 
time-asymptotic phase, the torque contributions become smaller and more localised. Transients remain important for the 
unrealistically slow turn on A = 1/2 (full width of 2 Gyr scaled to the Milky Way). The profiles for smaller, more realistic 
size bars are qualitatively similar are similar to those for the large bar (Panels d & e in Fig. QjJ ; ILR continues to dominate 
the total torque with corotation and (1, 0, 2) important at early times. For these smaller bars, the ILR changes sign after the 
main peak (T ~ 1) as the initial response laps the slowing bar before settling to the time-asymptotic LBK value. The torque 
is deposited as smaller radii and energy as expected but still over a broad region inside of the bar radius. Paper 2 presents 
estimates of the particle number required to simulate the resonant interaction in these examples. 



3.2 Decaying satellite 

3.2.1 Circular orbits 

Orbiting satellites couple to a live halo through resonances in much the same way as does a rotating bar. However, because 
the satellite size is generally much smaller than the halo scale, the total torque will have significant contributions from 
harmonic orders 2 larger than quadrupole. Just as in the case of the bar, the ( — 1,2,2) resonance can couple to very small 
radii although the satellite perturbation will have a different frequency spectrum and amplitude than the bar perturbation. 
Higher-order contributions have energies and response locations similar to the satellite. To investigate this situation, I assume 
that the satellite's orbit decays according the the local Chandrasekhar formula while remaining on a circular orbit. This, 
then determines the radius of the satellite orbit and the instantaneous orbital frequency, Qs- The low-order response is then 
evaluated using the theory from ij2| It is possible to compute the torque self-consistently by including a larger number of I 
terms (eg. 1141 . However, in this example, we only want to illustrate the contribution of the lower-order terms for comparison 
with the bar example and therefore the additional complexity and time to compute the self-consistent evolution is not justified. 
Because the satellite orbital radius changes in time, the torque takes the form of equation (1151 : the time-dependence in the 

2 The contribution scales as (lni) - 1 where I is the harmonic order (Weinberg 1986). 
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(e) a = 0.01, M b = 0.011, T a = 1/2, A = 1/4 



Figure 5. Relative contributions to the total torque L z = I z tl z as in Fig.[5]for a range of bar lengths (a), bar masses (Mg) and growth 
times (To and A). Only the top four contributors are shown individually. The total torque is the sum over all contributing terms for 
I = m = 2. 



© 2004 RAS, MNRAS 000.IT1IT51 



Time- dependent secular evolution 11 



0.22 
0.2 
0.18 
0.16 
0.14 
0.12 
0.1 
0.08 
0.06 
0.04 
0.02 




M s =0.1 
0.05 
0.025 
0.01 
0.0075 
0.005 
0.0025 



0.5 



1.5 



2.5 



3.5 



Figure 6. Change in radii with time due to orbital decay according to the local Chandrasekhar dynamical friction formula for nearly 
circular orbits for a variety of satellite masses M s (in units of the halo virial mass. 



rotation is not separable from the time-dependence of the changing Fourier coefficients in equation (JKJ. The torque is then 
straightforwardly computed for low-order I and m using equation I17|l although it is more computational intensive than the 
bar example due to the non-separable time dependence than the bar evolution case. 

The initial orbital radius is r = 0.2 (20% of the virial radius or approximately 60 kpc scaled to the Milky Way). Figure^] 
shows the radius as a function of time for a variety of satellite masses M 3 from 0.0025 to 0.1 virial mass units for interpreting 
the torque plots that follow. As in the case of the bar, the mass is slowly increased using f(t) (from H3.2I with To = 1/2 
and A = 1/4. The smallest mass satellite requires approximately 6 Gyr to decay from r = 0.2 and the most massive satellite 
decays in less that 100 Myr. 

Figure shows the torque, dL a /dt, as a function of time and orbital radius for the quadrupole I — m = 2 spherical 
harmonic contribution. The torque for all values of M 3 share some common trends: 1) the I = m = 2 torque is dominated by 
corotation (0:2:2). At early times, there is a moderate ILR transient (-1:2:2) which changes sign and oscillates at later times 
and the satellite position laps the transient wake, alternatively shifting the sign of the torque. We can remove the turn-on and 
extrinsic amplitude dependence by dividing dL s /dt by M% and multiplying the time by M s J* dt' f(t'). Each value of time, 
then, corresponds to the same satellite radius in each panel. The initial transient in dL z /dt between radii 0.12 and 0.2 (Panels 
b, d, and f of Fig.0 has a nearly self-similar form. The differences are due to intrinsic time dependence from transients. For 
example, as M s decreases, the ILR becomes more localised in time with a larger peak amplitude. 

Figure |H| illustrates the distribution of torque from the ILR as a function of time and energy and radius similar to Figure 
[3] Torque from the initial transient for the satellite orbiting at r ft* 0.15 ft* 2r 3 is deposited at radii between 0.005 and 0.02 
or between 7% and 30% of r s . Figure plots the ratio of the cumulative change in angular momentum in energy and time, 
AL(< E, < £)/£(< E), for the satellite decay shown in Figure^] As in the bar-halo interaction, the peak total torque occurs 
at much larger radii than the peak relative torque. Comparing to Figure [I] we see that the change in angular momentum is 
in the same region of phase space and larger than the 15%-20% fractional change found in the bar example. Since the halo 
evolves for the bar simulation, we anticipate evolution for the halo driven by this satellite. The N-body simulation is mor e 
difficult to perform in this case because of the disparate scales and will be the subject of a later paper l|Choi et alJbooA . 
Since the decay time is proportional to M s _1 and the torque is proportional to Ms, the total torque deposited in the transient 
phase is proportional to satellite mass . Therefore, as long as M s is sufficiently large that the orbital decay takes place in 
a galaxy age, the angular momentum transport from a few larger merger events may affect the halo profile as described in 
Paper 2. 
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Figure 7. Time-dependent torque for various resonances for a decaying circular satellite orbit with log A = 4 T = 1/2, A = 1/4, and 
initial radius Rq = 0.2. Only the top three contributing terms are shown individually. 
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Figure 8. ILR contribution to the torque as in Fig.E|but for the satellite decay with M s = 0.025. 
3.2.2 Eccentric orbits 

Eccentric orbits may be treated similarly using equation (1171) . Rather than enforce a circular orbit, one may solve the equations 
of motion for an orbit decaying according to the Chandrasekhar formula and use the solution, simultaneously, as the input 
to equation 11711 . In practice, it is convenient to solve the entire problem with a standard ODE solver. Unfortunately, there 
is a wide range of time scales in the problem: the decay time and the orbital time scales throughout phase space. In practice, 
then, accurate solutions require small time steps. 

The additional satellite orbital frequency complicates the spectrum and, because the orbital radius now oscillates as well as 
decays, complicates the representation of the interaction. However, qualitatively, the behaviour is similar to the circular case: 
low-order resonances such as the ILR provide strong positive torque on the halo at early times and then become oscillatory. 
Unlike the circular case, an eccentric orbit will deposit its orbital angular momentum at lower energies because the amplitude 
of the perturbation will be larger at smaller radii. This response of the halo may be a wake that influences the subsequent 
evolution of the disk and the halo. We will describe our detailed findings for the satellite-halo interaction for eccentric orbits 
in a later paper. 



4 SUMMARY 

This paper shows the lLvnden-Bell fc Kalnaisl l|l97l LBK) secular torque formula is not a quantitatively accurate description 
of secular evolution in galaxies. The LBK formula assumes that the growth time of a perturbation driving secular evolution 
is infinitely long compared the characteristic dynamical time. However, not only is the number of characteristic dynamical 
times in a galactic age modest, the growth of a spiral arms, a bar or the decay of a large dwarf satellite has evolutionary time 
scales of order 100 Myr to 1 Gyr, uncomfortably close to the equilibrium galaxy lifetime of 10 Gyr or smaller. We have seen 
that the finite time effects yield a different distribution of angular momentum in the halo and a different total torque than for 
a time-asymptotic system assumed by LBK. Similarl y, it is unlikely that time-dep e ndent effects will affect time- asymptotic 
perturbation theory predictions in planetary systems iGoldreich fc Tremainj dl979T t : IGoldreich fc Tremainel Jl98ut) where the 
secular evolution time scale is much larger than the dynamical time scale. 

In |J5| I presented generalised the LBK formula for perturbations of finite duration. The new formula shares many of the 
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Figure 9. Ratio of cumulative ILR angular momentum transfer to total angular momentum as in Fig.|l]but for the satellite decay with 
M a = 0.025. 



features of the LBK formula; the main difference is that the Dirac delta function 8(1- fl — mQ p ) is replaced with a time integral 
(see ea. 1151 . Over astronomically- relevant time scales, the finite-time generalisation shows that transients from the formation 
history of the system play an important role in the overall dynamical evolution. For example, Figures and [5] illustrate the 
total torque and differential contributions of the principal resonances during slow down of a bar in dark matter halo. The 
correct result is the net contribution for approximately 8 resonances (not all shown in Fig. |5J and therefore these resonances 
must be accurately represented in an N-body simulation to obtain the correct result (see Paper 2). 

The most important general finding is that the history of galaxy evolution can not be ignored in understanding and 
predicting a particular galaxy's evolution. These dynamics are illustrated in two cases: the bar-halo interaction and the 
satellite-halo interaction. The details of the bar-halo interaction have been rece ntly described in Paper 2 and explored in 
a self-consistent simulation in Paper 3 l)Kellv Hollev-Bockelmann fc Katd 12004^1 . Both papers demonstrate that bar-halo 
coupling through the ILR agrees with the predictions presented here. Both the LBK formula and the time-dependent formula 
predict that the ILR resonance dominates the overall angular momentum transfer. However, the time-dependent evolution 
spreads the resonant interaction over a broader range of lower energies that are more populated in phase space and therefore 
have a larger overall evolutionary consequence. The time-asymptotic ILR is located at smaller energies and radii and would 
be inaccessible to most N-body simulations. The time-asymptotic approximation is especially poor for the ILR resonance in 
a cuspy halo where the time-dependent spread in frequency corresponds to a large region in radius and energy. This is strong 
encouragement to attempt an understanding of the dynamical interactions prior to simulation. One expects a similar halo 
response in the satellite-halo interaction because the quadrupole perturbation from a rotating bar shares much in common 
with the quadrupole perturbation from an orbiting satellite. Indeed, we show that this is the case. However, because the 
symmetries, the radial profile, and the distribution of frequencies for the two perturbations are different, the ILR no longer 
dominates in the satellite interaction. Nonethel ess, the ILR does play a role during the early stages of orbital evolution and 
this has been verified in N-body simulation (see lChoi et al.ll2004l) . The aggregate effect of this mechanism from a population 
of substructure will be the subject of a later investigation. 

The wakes excited by interactions between satellites, disks and dark matter halos give rise to a variety of observable 
consequences. For example, large-scale transport of angular momentum throughout the disk and to the halo may will tend to 
decrease the disk scale length. The response of the galaxy halo to the surrounding group environment can propagate features to 
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the inner galaxy, exciting disk waves and subsequent star formation. In addition, the interaction between the satellite and disk 
will drive bending modes and the energy in these modes will heat the disk. Detailed simulations for these processes of disparate 
scales are difficult. Time-asymptotic secular perturbation theory has often been used to estimate evolution and compare with 
simulation. However, we have found that the difference in amplitude and significance between the time-dependent and time- 
asymptotic responses for the bar-halo interaction (Fig.|2J is dramatic. If these differences are representative of the difference 
in these other processes, our current picture of secular evolution must be revisited by considering the time-dependent evolution 
of these "waves" in the context of a galaxy's environmental history. 
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